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j  Abstract 

A  numerical  method  to  solve  the  perfect  gas  Navier- Stokes  equations  for  hypersonic 
flows  past  three-dimensional  blunt  bodies  has  been  developed.  The  numerical  method 
uses  flux-splitting  and  shock-fitting  with  an  implicit  Gauss-Seidel  line-relaxation  pro¬ 
cedure  to  accelerate  convergence.  The  technique  has  been  used  to  solve  the  flow  field 
over  a  spherically  blunted  biconic  and  the  X24C-10D  hypersonic  research  vehicle.  The 
method  has  been  shown  to  reduce  the  number  of  iterations  required  to  achieve  conver¬ 
gence  of  a  typical  problem  by  a  factor  of  about  one  hundred  over  an  explicit  method. 
The  scheme  also  shows  a  potential  advantage  over  approximately  factored  implicit 
methods.  The  key  advantage  of  this  technique  is  that  the  low  number  of  iterations 
required  for  convergence  does  not  increase  as  mesh  resolution  is  refined.  t\  ■ 
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A,B,C  =  Jacobians  of  F,  G  and  H 
c  =  local  speed  of  sound 

cp  =  specific  heat  at  constant  pressure 

cv  —  specific  heat  at  constant  volume 

e  =  total  energy  per  unit  mass 

ei  =  internal  energy  per  unit  mass,  e*  =  evT 

F,G,H  =  flux  vectors  of  the  conserved  flow  variables 
F  =  vector  form  of  the  flux  components 

i,j,k  =  mesh  coordinates  in  the  £,q,  ?  directions 
I,J,K  =  largest  values  of  the  mesh  coordinates 
k,  kt  =  molecular  and  eddy  conductivity  of  the  fluid 
M  =  Mach  number 

n  =  time  level  of  computation 

p  =  pressure 

Pr,  Prt  =  Prandtl  number  and  turbulent  Prandtl  number 
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=  heat  flow  vector 
=  nose  radius 
=  surface  normal  vector 
=  x  component  of  i  surface  normal  vector,  etc. 

=  Stanton  number,  St  =  |n  •  ArVr|/poo«oo(cproo  +  -  cpTw*\\) 

=  time 

=  temperature 

=  x,  y  and  z  components  of  velocity 
=  vector  of  conserved  flow  variables 
=  cell  volume 

=  vector  of  the  primitive  flow  variables 
=  coordinate  directions 
=  law  of  the  wall  variable 
=  angle  of  attack 

=  implicit  change  in  U,  6Un  =  Un+1  -  Un 
=  explicit  change  in  U,  defined  in  Equation  (7) 

=  ratio  of  specific  heats 
=  transformed  coordinates 
=  molecular  and  eddy  viscosity  coefficients 
=  bulk  viscosity  coefficient 
=  fluid  density 
=  shear  stress  tensor 
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Introduction 


The  numerical  solution  of  the  three-dimensional  Navier-Stokes  equations  is  becoming 
commonplace  and  large  complicated  problems  have  been  solved.  A  notable  example 
is  the  solution  of  the  flow  field  over  an  entire  hypersonic  aircraft  performed  by  Shang 
and  Scherr1  in  1985.  The  numerical  solution  of  the  Navier-Stokes  equations  is  a  com¬ 
putationally  intensive  task  requiring  the  full  resources  of  a  supercomputer.  A  typical 
problem  may  require  on  the  order  of  104  iterations  to  achieve  a  steady-state  result  when 
using  an  explicit  algorithm.  This  is  a  result  of  the  strict  time  step  limitation  imposed 
on  the  solution  because  of  the  fine  mesh  spacing  needed  near  surfaces  to  resolve  viscous 
and  heat  transfer  effects.  The  Navier-Stokes  equations  also  require  a  large  number 
of  calculations  per  time  step  to  compute  all  of  the  terms  in  the  equations.  Thus  the 
combination  of  these  two  factors  leads  to  the  use  of  much  computer  time  to  achieve 
a  steady-state  solution.  A  numerical  method  that  reduces  the  number  of  iterations 
required  to  reach  a  steady-state  or  reduces  the  number  of  computations  required  per 


iteration  would  enhance  the  ability  to  study  flow  physics  using  numerical  solutions  or 
allow  more  test  cases  to  be  run  for  a  given  amount  of  computer  time. 

In  1985  MacCormack2  presented  a  numerical  algorithm  for  the  solution  of  the  two- 
dimensional  Navier-Stokes  equations  that  demonstrates  a  large  improvement  in  the 
number  of  iterations  required  for  a  steady-state  result.  The  method  was  tested  on 
several  problems  and  a  solution  in  less  than  twenty  iterations  was  obtained  for  a  two- 
dimensional  transonic  nozzle.  The  three-dimensional  analog  of  this  technique  has  been 
developed  and  tested;  its  characteristics  and  some  test  results  from  this  method  are 
presented  here.  The  numerical  method  demonstrates  several  positive  qualities,  among 
which  it: 

•  is  fully  conservative  to  prevent  any  numerical  loss  of  mass,  momentum  or  energy. 

•  is  fully  implicit,  allowing  large  time  steps  to  be  taken. 

•  uses  Gauss-Seidel  line-relaxation  to  accelerate  convergence. 

•  uses  flux  dependent  differencing  to  maintain  stability  without  having  to  introduce 
any  smoothing  terms  into  the  difference  equations. 

•  is  second-order  accurate  in  all  spatial  directions. 

•  uses  shock-fitting  to  accurately  capture  the  very  strong  gradients  in  the  flow  vari¬ 
ables  that  occur  at  the  bow  shock  wave. 

However,  the  numerical  method  does  have  one  main  drawback,  it  requires  many  com¬ 
putations  per  time  step,  so  that  a  large  reduction  in  iterations  converts  to  a  smaller 
reduction  in  computer  time.  The  results  presented  were  obtained  using  a  program  that 
is  not  vectorized;  a  further  reduction  in  computer  time  would  be  realized  if  vectorization 
was  performed. 

The  numerical  method  has  been  tested  on  two  configurations  and  the  results  compare 
well  with  previous  experimental  and  computational  work.  The  first  body  tested  was  a 
spherically  blunted  cone  travelling  at  a  Mach  number  of  7.97  and  a  Reynolds  number  of 
1.2  x  10g/m;  angles  of  attack  of  between  0.0°  and  7.0°  were  used.  Secondly,  the  flow  over 
the  X24C-10D  hypersonic  research  vehicle  at  a  Mach  number  of  5.95,  a  characteristic 
Reynolds  number  of  16.4  x  106/m,  and  an  angle  of  attack  of  6.0°  was  computed.  This 
is  a  duplication  of  the  explicit  efforts  of  Shang  and  Scherr,  but  using  implicit  Gauss- 
Seidel  line-relaxation  auad  shock-fitting.  The  use  of  the  implicit  method  allowed  a  better 
resolution  of  the  near  wall  viscous  and  heat  transfer  effects.  The  number  of  time  steps 
required  to  achieve  a  steady-state  solution  for  these  problems  is  about  a  factor  of  one 
hundred  less. 

The  Navier-Stokes  Equations 

The  Navier-Stokes  equations  may  be  written  as  the  integration  over  the  surface  of 
an  arbitrary  volume 
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where  0  is  the  vector  of  conserved  variables  averaged  over  the  volume,  F  is  the  vector 
form  of  the  fluxes  of  these  conserved  quantities,  F  =  Fit  +  GTy  +  Hi,,  and  dS  is  the 


differential  surface  vector  of  the  volume.  F,  G,  and  H  are  the  flux  vectors  expressed  in 
Cartesian  coordinates  given  below, 
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The  equation  set  is  closed  by  using  the  perfect  gas  law  to  relate  temperature  to 
pressure  and  density,  the  Sutherland  viscosity  formula  for  air3  to  determine  viscosity 
as  a  function  of  temperature,  the  Stokes  relation  to  give  the  bulk  viscosity,  and  the 
assumption  of  constant  Prandtl  number  to  determine  the  conductivity  of  the  gas.  The 
two  turbulent  transport  coefficients,  nt  and  fc*,  are  determined  using  the  Baldwin- 
Lomax  turbulence  model4. 
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where  the  subscript  r  represents  a  reference  temperature  and  Ti  is  taken  to  be  110°K. 
These  equations  were  used  for  the  present  calculation  in  which  temperatures  where  well 
below  those  where  dissociation  or  ionization  of  the  gas  may  take  place. 

The  set  of  differential  equations  is  fully  described  by  imposing  appropriate  initial 
and  boundary  conditions.  The  initial  conditions  are  not  critical  because  the  required 
solution  is  steady-state,  therefore  an  arbitrary  initial  condition  amounting  to  a  rough 


guess  of  the  initial  now  field  and  satisfying  the  constitutive  relations  can  be  imposed. 
The  boundary  conditions  are  determined  by  the  physical  nature  of  each  boundary.  At 
the  surface  of  the  body  a  no-slip  condition  is  used  with  either  an  adiabatic  wall  or  a 
fixed  wall  temperature.  At  the  boundary  outside  of  the  bow  shock  wave,  which  is  the 
free-stream,  the  flow  conditions  are  fixed  since  they  are  uninfluenced  by  the  interior 
flow  field.  In  the  streamwise  direction  the  only  relevant  boundary  is  that  at  the  end  of 
the  body;  this  is  set  by  assuming  that  =  0,  which  is  physically  consistent  because 
the  flow  is  either  supersonic  or  within  the  boundary  layer  at  this  point.  The  boundaries 
in  the  meridional  direction  are  set  using  bilateral  symmetry  considerations. 

The  Numerical  Method 

The  numerical  procedure  used  is  a  generalization  to  three  dimensions  of  the  implicit 
Gauss-Seidel  line- relaxation  technique  developed  by  MacCormack  for  two  dimensions 
and  is  derived  from  the  volume  formulation  of  the  Navier-Stokes  equations.  The  first- 
order  difference  equations  are  derived  here,  although  the  generalization  to  second-order 
spatial  accuracy  follows  directly. 

Consider  a  three-dimensional  mesh  made  up  of  finite  volumes  with  the  flow  variables 
stored  at  the  centroids  of  each  volume,  and  the  corners  of  each  volume  stored  as  the 
i,  y,  z  coordinates  of  the  grid  points.  We  can  think  of  the  discrete  form  of  the  volume 
formulation  of  the  Navier-Stokes  equations  as  a  summation  of  the  fluxes  entering  a 
finite  volume: 
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l2> 

V 

We  now  adopt  a  sign  convention  on  S  so  that  a  component  of  S  is  positive  when  directed 
in  increasing  £,  r],  or  f.  We  also  break  F  into  three  parts  as 

F  =  /+  +  F-  +  /. 

where  F+  represents  the  flux  induced  by  the  inviscid  terms  in  the  equations  moving  in 
the  positive  £,  rj,  or  f  direction  and  F_  is  the  negatively  moving  flux.  F  is  the  flux 
produced  by  the  viscous  stresses  and  thermal  effects.  Discretizing  Equation  (2)  we 
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where  S,+  i  denotes  the  surface  vector  of  the  surface  between  point  i,j,k  and  point 
i  +  \,j,k,  with  the  other  surface  vectors  defined  similarly.  (Note  that  because  of  the 
definition  of  S  the  above  summation  is  made  up  of  positively  and  negatively  signed 
components.)  We  let 
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with  the  other  components  of  the  fluxes  defined  similarly.  Linearizing  the  flux  vectors 
in  time  gives 

F'T„\  =  Pin*  +  A'lxjkSU:]k  +  0(At2), 

= G>n+iJk + B'%ksu:]k + o^), 

H%\  =  H'%k  +  C"ll}kbu:jk  +  0(At2). 

Where  A',  B',  and  C'  are  the  Jacobians  of  the  flux  vectors  of  F' ,  G',  and  H'  with  respect 
to  U ,  and  6Un  =  Un+l  -  Un.  Here  we  are  using  first-order  accuracy  in  time  because  we 
are  interested  in  the  steady-state  solution  and  temporal  accuracy  is  not  required.  The 
method  may  be  made  second-order  accurate  in  time  with  some  effort.  The  finite  volume 
Navier-Stokes  equations  may  now  be  written  more  succinctly  using  difference  operators 
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We  let 


F'ZV  =  Kk  +  *K- 


with  similar  expressions  for  the  other  two  viscous  fluxes.  We  apply  the  thin-layer 
assumption  to  the  derivatives  of  the  implicit  viscous  terms.  This  implies  that  the 
differences  with  respect  to  r\  are  much  larger  than  those  with  respect  to  £  or 
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Thus  we  need  only  find  an  expression  for  6G’. 

We  define 
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to  obtain 


G'  =  FS'JZ  +  GS'jy  +  HS 


The  rotated  viscous  flux  G'  is  made  up  of  stresses  that  contain  derivatives  with  respect 
to  the  Cartesian  variables  z,y,z  which  must  be  changed  to  the  general  variables 
However,  if  we  apply  the  thin  layer  assumption,  the  derivatives  with  respect  to  £  and 
<;  may  be  neglected.  Therefore  we  have 
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The  derivatives  of  the  mesh  variables  with  respect  to  the  Cartesian  variables  are  de¬ 
termined  using  the  relations  derived  from  the  coordinate  transformation  matrices 
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Finally,  using  these  changes  of  variable,  we  can  rewrite  6G'  as 

SC'-  3  fj5Vl 

6G  - - —  Mr,  —  (6V), 

where 

Mr,  = 


o  5/1(2m+A)+m(5/v  +  S/,) 
0  Sj*SjV(/i+A) 

0  S]sS]t(ji  +  \) 


SjTSjy(i*+*) 
S;v(2A+A)+M(S/t  +  S/j 


Sjv5„(m+A) 
S/i(2M  +  A)  +  M(SJ2I+S;2y) 


"*5,2  =  u(SjX(2ii  +  A)  +  (i(Sfy  +  S2j)  +  vSJZS}y(p,  +  A)  +  +  A), 

m5  s  =  uSJXS]y(ti+  A)  -(-  v[Sjy(2ii  +  A)  +  t(S*x  +  S2j)  +  wS]vS]x(m  +  A), 
ms  4  =  uS]ZS]x(fi  +  A)  +  vS)VSJX{ii  +  A)  +  w(S}x(2 m  +  A)  +  ^(S2.,  +  5;2V))- 

"*5.s=  -(5/,+5/v+S*,). 
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Here  the  bars  on  the  velocities  represent  the  average  value  of  that  velocity  between 
the  two  volumes  in  question.  We  can  change  variables  from  6V  to  6U  by  defining  the 
Jacobian  N  to  be  N  =  so  that  <5V  =  N6U.  Thus  we  finally  have 

Combining  Equations  (4)  and  (6),  we  have  the  final  form  of  the  finite  volume  difference 
equation  which  approximates  the  Navier-Stokes  equations.  The  equation  is  given  below 
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with  the  quantities  that  are  known  at  time  n  on  the  right  hand  side,  and  those  that 
must  be  solved  implicitly  at  time  n+1  on  the  left. 
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(7) 


This  completes  the  derivation  of  the  difference  equation  used  to  solve  the  Navier-Stokes 
equations. 

The  implementation  of  the  boundary  condition  for  the  explicit  side  of  the  difference 
equation  is  straightforward.  A Un  is  computed  using  the  prescribed  boundary  conditions 
mentioned  above.  However,  the  implicit  boundary  conditions  are  more  complicated. 
Let  us  write  Equation  (7)  in  terms  of  5x5  block  matrices 


A6U:]k  +  B6U”  +  lk  +  C6U»_lk  +  D6U?+l]k  +  E6U?_l)k  +  F6U?,k. M  +  GSU^.j  =  A  t/,%,  (8) 

where  the  matrices  with  tildes  are  defined  implicitly  from  Equation  (7).  We  may  write 
Equation  (8)  in  matrix  form  at  a  general  i  and  k  location  as 
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The  top  and  bottom  equations  are  used  to  set  the  freestream  and  the  wall  boundaries 
respectively.  The  freestream  boundary  is  easy  to  set  because  it  is  a  supersonic  inflow 
so  that  the  flow  field  remains  invariant.  Thus  we  have  6UtJk  =  o,  which  implies  that 
Ej  =  o,  and  therefore  Bj~  j  is  redundant.  For  the  wall  boundary  condition,  we  impose 
that  there  is  no-slip  at  the  wall,  the  normal  derivative  of  pressure  at  the  wall  is  zero 
(from  boundary  layer  theory),  and  that  a  wail  temperature  condition  holds  depending 
on  whether  the  wall  has  a  fixed  temperature  or  is  adiabatic.  These  conditions  at  the 
wall  fix  the  flow  field  at  the  wall  uniquely  and  imply  that  5u',  =  -6u2,  bv\  =  - 6v'2 , 
Sw\  =  -6ui'2,  and  6pi  =  Spi.  For  the  internal  energy,  we  impose 
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adiabatic  wall; 
fixed  wall  temperature. 


Here  we  are  using  a  grid  where  a  dummy  point  within  the  body  is  used  (subscript 
j  =  1)  to  set  the  boundary  at  the  body  surface.  The  primes  on  the  velocities  represent 
rotated  velocities  given  by 

ti1 2 3  =  uStx  +  vSty  +  u/Stz, 
v'  =  u  SJZ  +  vS]y  +  wSjz, 
w'  ~  u Skx  +  vSk y  +  wSkz. 

We  can  write  the  wall  boundary  condition  in  matrix  form  as  Ti 6UX  =  T26U2  which 
implies  that  E2  =  T^Tj.  Therefore  the  next  to  last  equation  in  (9)  may  be  written  as 


B26Utsit  +  -t- 


We  define 


A'2  =  ( A2  +  C2E1 ). 

Then  by  contracting  the  upper  and  lower  equations,  we  rewrite  (9)  as 


Aj-i  Cj-i 


6Uij- 1 


Bt  A'2 
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A  similar  treatment  of  the  boundaries  in  the  1  and  k  directions  for  the  implicit  side 
of  the  Equation  (7)  is  required.  The  procedure  generalizes  from  the  above  discussion 
directly. 

This  completes  the  derivation  of  the  technique  to  solve  the  Navier-Stokes  equations 
using  Gauss-Seidel  line-relaxation.  The  procedure  to  solve  these  difference  equations 
is  as  follows: 

1.  Solve  the  right  hand  side  of  Equation  (7),  for  the  whole  flow  field  using 

data  at  the  current  time  level  n. 

2.  Use  A as  a  driving  term  to  solve  the  left  hand  side  of  the  equation,  using 
the  latest  available  data  (i.e.  using  a  Gauss-Seidel  technique).  This  is  done  by 
solving  a  block  tridiagonal  matrix  equation  in  6U*k  along  a  particular  j  line  in 
the  flow  field.  The  solution  line  is  swept  through  the  flow  field  from  the  farthest 
downstream  t  plane  to  the  upstream  i  plane,  then  the  reverse  is  performed.  The  j 
solution  line  is  swept  through  all  k  locations  during  each  streamwise  pass.  Thus, 
with  a  large  enough  time  step,  each  point  in  the  flow  field  could  conceivably  be 
influenced  by  every  other  point. 

3.  Update  the  solution  using  and  continue  the  process  until 

convergence  is  achieved. 

Implementing  second-order  spatial  differencing  is  fairly  straightforward  and  is  accom¬ 
plished  by  replacing  for  example  F  ^  =  F*,  (first-order),  with  F  x  =  |F+,  -  |F.,_  1 


(second-order).  In  the  j  direction  on  the  implicit  side,  the  data  from  the  most  re¬ 
cent  time  step  is  used  to  compute  the  fluxes  originating  from  two  cells  away  from  the 
cell  in  question.  This  maintains  the  block  tridiagonal  structure  of  the  matrix  solution 
without  sacrificing  second-order  accuracy.  Although  the  fluxes  from  two  cells  away 
are  somewhat  lagged  in  time,  this  becomes  unimportant  as  the  solution  converges  to  a 
steady-state. 


The  Shock  Fitting  Routine 

The  technique  of  shock-fitting  was  originally  developed  by  Moretti5,  however  the 
method  used  for  this  solution  is  somewhat  different  and,  as  a  result,  is  presented  here. 
The  bow  shock  wave  that  envelops  a  hypersonic  body  is  represented  as  a  surface  of 
discontinuity  which  is  moved  through  the  flow  field  in  a  physically  consistent  manner. 
Computationally,  the  shock  wave  is  located  at  a  constant  j  mesh  surface  and  acts  as 
a  boundary  in  the  solution.  The  initial  position  of  the  shock  wave  is  guessed  at  the 
beginning  of  the  calculation,  and  then  the  shock  wave  is  moved  through  the  mesh 
after  each  time  step  until  the  steady-state  position  is  reached.  The  procedure  for  each 
iteration  is  to  solve  the  discretized  Navier-Stokes  equations  for  the  entire  flow  field. 
Then  the  shock  wave  velocity  is  computed  at  the  mesh  surface  that  represents  the 
s‘  ^ck  wave.  The  mesh  is  then  reconfigured  by  moving  the  shock  wave  surface  inward 
or  outward,  and  all  of  the  other  interior  grid  points  are  adjusted  in  a  consistent  manner. 

The  shock  wave  velocity  at  each  point  is  found  by  solving  two  of  the  Rankine- 
Hugoniot  shock  jump  relations  for  an  oblique  shock  wave  and  one  characteristic  re¬ 
lation  for  a  backward-moving  characteristic.  We  denote  freestream  conditions  with  a 
subscript  0,  the  conditions  at  a  location  immediately  behind  the  shock  with  a  subscript 
shk,  and  the  conditions  at  the  centroid  of  the  first  volume  behind  the  shock  surface 
with  a  subscript  1.  This  nomenclature  is  used  because  the  shock  surface  is  located 
between  the  centroids  of  two  finite  volumes,  i.e.  there  is  a  finite  distance  between  the 
shock  surface  and  the  first  downstream  point  where  data  exists.  The  characteristic 
relation  is  used  to  tie  the  conditions  at  these  two  locations  together.  The  magnitude  of 
the  shock  wave  velocity  is  denoted  as  u;shk  and  has  direction  parallel  to  the  freestream 
velocity. 

The  relevant  shock  jump  equations3  for  a  shock  moving  with  speed  u'shk  are 
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The  relation  for  a  backward  running  characteristic  is 

UnSHK  -  Unl  =  — (PSHK  ~  Pi) 
PC 
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The  element  of  the  shock  surface  under  consideration  is  moved  a  distance  wshkAI  in  the 
direction  parallel  to  u0.  The  point  where  the  plane  formed  by  the  new  shock  surface 
intersects  the  surrounding  mesh  lines  is  found,  the  interior  of  the  mesh  is  adjusted, 
and  the  change  in  volume  of  each  volume  is  saved  to  account  for  the  unsteady  mesh 
transformation  terms  in  the  calculation  of  the  surface  fluxes. 

It  should  be  noted  that  the  shock-fitting  technique  has  the  advantage  of  satisfying 
the  Rankine-Hugoniot  shock  relations  exactly  at  a  single  surface  in  the  computation, 
rather  than  smearing  the  shock  wave  across  several  grid  surfaces  as  in  a  shock  capturing 
method.  Thus  fewer  mesh  surfaces  are  needed  in  the  vicinity  of  the  bow  shock  wave. 
Any  imbedded  shock  waves  in  the  flow  field  such  as  the  canopy  and  strake  induced 
shocks  are  captured  by  the  numerical  method. 


Results  of  Test  Cases 

Biconic  Results 

The  numerical  method  was  tested  on  two  different  configurations.  The  first  was  a 
spherically  capped  biconic  body  travelling  through  air  at  a  Mach  number  of  7.97  with  a 
characteristic  Reynolds  number  of  1.2  x  108/m.  This  configuration  is  of  interest  because 
it  has  been  studied  numerically  and  experimentally2,6,7,8.  It  consists  of  a  spherical  cap 
of  radius  rn  joined  to  a  cone  of  half  angle  10.5°,  which  joins  a  second  cone  of  half  angle 
of  7.0°  at  a  point  21.5  nose  radii  from  the  tip.  Bilateral  symmetry  was  assumed  and  a 
body-fitted  mesh  was  generated  algebraically  using  exponential  stretching  away  from 
the  body  surface  and  along  the  cone.  The  initial  bow  shock  wave  location  was  guessed 
to  be  a  paraboloid  of  revolution,  and  then  the  shock  was  allowed  to  move  during  the 
solution  as  described  above.  A  cross-section  of  the  initial  mesh  through  the  bilateral 
plane  of  symmetry  is  shown  in  Figure  1. 

The  solution  was  performed  in  two  sections,  first  a  converged  result  for  the  nose 
section  (about  4rn  long)  was  obtained.  Then  this  solution  was  used  as  an  inflow 
boundary  for  the  afterbody  section  of  the  calculation.  This  is  physically  correct  because 
the  flow  is  purely  supersonic  or  within  the  boundary  layer  at  the  plane  where  the 
solution  was  broken.  This  solution  technique  is  more  efficient  than  performing  the 
entire  calculation  at  once  because  the  flow  field  around  the  afterbody  is  dependent  on 
that  around  the  nose  section.  Therefore,  until  the  solution  for  the  nose  is  converged,  the 
afterbody  solution  will  not  progress  to  the  correct  solution  and  cannot  converge.  The 
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converged  solution  was  obtained  in  120  iterations  for  the  nose  section  on  a  mesh  with 
30  points  in  the  streamwise  direction,  20  in  the  direction  normal  to  the  body,  and  18  in 
the  meridional  direction.  The  results  for  the  afterbody  were  computed  in  180  iterations 
on  a  18x20x18  mesh.  A  maximum  CFL  number  of  1500  was  used.  The  converged  grid 
shape  is  shown  in  Figure  2  with  the  shock  wave  position  indicated.  Clearly  the  grid 
has  been  deformed  greatly  and  has  moved  a  significant  physical  distance.  A  plot  of 
the  surface  pressure  distribution  on  the  afterbody  compared  to  experimental  evidence 
from  Coutler  and  Carver8  is  given  in  Figure  3.  The  results  show  good  agreement  with 
experiment,  with  a  maximum  deviation  of  about  3%.  A  plot  showing  velocity  profiles 
at  two  x/rn  locations  is  presented  as  Figure  4;  it  shows  the  well  defined  boundary  layer 
which  has  been  resolved  by  the  method  (the  y+  at  the  first  point  from  the  wall  ranged 
from  0.3  to  1.2  in  this  case).  Figure  5  is  a  plot  of  the  temperature  profile  at  the  same 
two  x/r„  locations.  It  should  be  noted  that  the  solution  is  for  a  fixed  wall  temperature 
of  311.1  °K.  The  figure  shows  a  very  good  resolution  of  the  very  large  temperature 
gradient  near  the  wall  and  the  inflection  point  in  the  temperature  distribution  caused 
by  the  interaction  of  the  cooled  wall  with  the  high  temperature  near-wall  stream.  The 
shock  standoff  distance  computed  with  the  current  method  agrees  with  the  results  of 
Scherr  and  Shang7  to  three  significant  figures. 

A  second  solution  of  the  zero  angle  of  attack  cone  was  performed  with  the  same 
number  of  points  in  the  grid,  but  with  a  near-wall  spacing  3.0  times  finer  than  the 
previous  result  (the  y+  for  the  first  point  from  the  wall  was  reduced  by  a  factor  of 
2.9).  The  boundary  layer  and  thermal  layer  are  essentially  identical  for  each  case.  The 
number  of  time  steps  required  to  reach  a  steady-state  remained  the  same. 

The  flow  field  over  the  same  body  at  a  7.0°  angle  of  attack  with  the  same  flow 
conditions  was  also  computed.  The  solution  was  obtained  in  135  iterations  for  the  nose 
section  and  200  iterations  for  the  afterbody.  This  is  about  a  10%  increase  over  the 
zero  angle  of  attack  case.  The  final  grid  is  given  in  Figure  6  and  shows  that  the  bow 
shock  wave  is  swept  more  closely  to  the  body  on  the  windward  side  than  on  the  leeward 
side.  The  asymmetry  is  also  evident  in  the  surface  pressure  plot  for  the  windward  and 
leeward  sides  of  the  body,  which  is  given  in  Figure  7. 


X24C  Hypersonic  Vehicle  Results 


The  flow  field  over  the  X24C-10D  hypersonic  vehicle  was  computed  for  free-stream 
conditions  of  M =  5.95,  an  angle  of  attack  of  6.0°,  and  a  characteristic  Reynolds  num¬ 
ber  of  16.4  x  106/m.  The  X24C  was  a  proposed  configuration  for  a  second  generation 
space  shuttle  which  has  a  spherically  blunted  nose,  a  canopy,  and  wings  and  strakes 
in  the  aft  section  of  the  aircraft.  This  body  is  a  difficult  test  case  for  the  numerical 
method  with  several  complicated  physical  phenomena  occurring  in  the  flow  field.  The 
grid  used  was  provided  by  Shang  and  Scherr1  and  is  comprised  of  60  x  18  x  28  points 
for  the  nose  and  midsection  of  the  body,  and  58x18x55  points  for  the  afterbody,  or  a 
total  of  87,660  grid  points.  The  surface  geometry  of  the  X24C  is  presented  in  Figure  8. 
In  this  case  the  solution  was  broken  into  four  sections  to  improve  the  efficiency  of  the 
solution  and  also  to  allow  each  section  to  fit  in  core  memory  (i.e.  less  than  two  million 
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words  available  on  a  Cray  X-MP  48). 

The  results  show  good  agreement  with  the  results  of  Shang  and  Scherr1  and  the  ex¬ 
perimental  results9,10,11.  The  numerical  method  accurately  predicts  the  surface  pres¬ 
sure  on  the  vehicle  except  in  the  tail  region  on  the  leeward  side,  as  seen  in  Figure  9.  The 
poor  resolution  of  the  tail  surface  pressure  is  probably  caused  by  the  coarse  grid  near 
the  tail  which  is  shown  in  Figure  10.  Shang  and  Scherr  used  twice  as  many  grid  points 
in  the  meridional  direction  and  were  able  to  predict  the  pressure  more  accurately,  but 
a  finer  grid  could  not  be  used  in  this  study  because  of  the  memory  limitation.  Figure 
10  also  demonstrates  that  the  numerical  technique  does  not  require  many  mesh  points 
in  the  flow  field  where  the  gradients  are  relatively  small.  Also  the  shock-fitting  method 
will  resolve  the  strong  bow  shock  wave  with  a  coarse  mesh  in  the  vicinity  of  the  shock 
wave.  This  results  in  an  efficient  use  of  grid  points  in  the  normal  direction.  Figure  11 
shows  the  Mach  number  distribution  on  the  plane  of  symmetry,  including  the  shock 
wave  and  recompression  of  the  flow  caused  by  the  canopy. 

The  next  five  figures  are  for  data  at  the  x/rn  =  108  plane  where  much  of  the  experi¬ 
mental  data  was  taken9,10,1 1 .  The  first  plot  is  a  pitot  pressure  survey  with  experimental 
points10  indicating  the  location  of  the  enveloping  shock  wave.  The  agreement  with  ex¬ 
periment  is  very  good,  and  the  results  also  compare  well  with  those  of  Shang  and 
Scherr1.  Figure  14  is  a  plot  of  the  grid  at  x/rn  =  108,  which  shows  the  location  of  the 
fitted  bow  shock  surface  in  relation  to  the  experimental  results10.  The  agreement  is 
very  good,  indicating  that  the  shock-fitting  routine  produces  an  accurate  shock  location 
even  where  the  shock  is  oblique.  Figure  14  shows  the  surface  pressure  distribution  ver¬ 
sus  the  normalized  arc  length,  f/fMAX*  at  this  x/rn  location.  (f/fMAx  =  o  is  the  leeward 
point  and  r/fMAx  =  1  is  the  windward  point.)  The  numerical  results  agree  very  well 
with  experiment9.  The  peak  in  pressure  at  c/jmax  =  0-58  is  caused  by  the  protrusion 
of  the  leading  edge  of  the  body  into  the  flow.  The  sudden  pressure  drop  to  the  lee 
of  this  point  is  a  result  of  the  expansion  of  the  flow  as  it  accelerates  past  the  leading 
edge.  Figure  15  is  a  plot  of  log  St  (a  non-dimensional  measure  of  the  heat  transfer) 
which  shows  at  present  a  less  accurate  prediction  of  experimental  results.  The  heat 
transfer  is  resolved  well  on  the  windward  surface,  but  the  results  show  poor  agreement 
on  the  leeward  surface.  Figure  16  is  a  streamline  trace  on  the  leeward  surface  of  the 
body.  It  shows  that  there  is  a  large  cross-flow  separation  bubble  on  the  leeward  surface 
at  this  x/rn  location.  It  is  likely  that  the  grid  spacing,  though  able  to  resolve  this 
structure  in  the  flow  and  the  surface  pressure  distribution,  is  inadequate  to  capture  the 
temperature  gradients  in  this  region.  The  large  dip  in  the  Stanton  number  that  occurs 
at  f/fMAX  =  0.06  is  located  at  the  reattachment  point  of  the  vortical  structure  on  the 
leeward  surface.  The  heat  transfer  at  this  point  is  very  small  as  would  be  expected  due 
to  the  reduction  in  the  velocity  near  the  surface. 

Figures  17  and  18  are  streamline  traces  on  the  leeward  surface  of  the  body  at  the 
constant  x/rn  =  127.5  plane  where  the  central  fin  and  the  strake  may  be  seen.  The 
vortical  structure  that  was  present  at  x/rn  =  108  has  become  more  pronounced  and 
the  separation  region  covers  a  larger  portion  of  the  body.  A  secondary  vortex  is  evident 


between  the  main  vortex  and  the  wail.  Figure  19  is  an  enlargement  of  this  structure 
and  shows  the  level  of  complexity  of  the  flow  field  that  has  been  resolved.  Figure  18  is 
a  streamline  trace  at  x/rn  =  140  where  the  fin  and  strake  are  at  their  maximum  size. 
The  vortical  structure  may  be  seen,  but  its  strength  has  diminished  because  the  strake 
has  reduced  the  cross-flow  which  drives  the  vortex. 

The  number  of  iterations  required  to  reach  a  steady-state  solution  varied  as  presented 
in  the  table  below. 


t  Range 

Iterations 

1-  24 

190 

24-60 

190 

60-84 

120 

84  -  118 

250 

Table  1.  Iteration  count  for  each  solution  section. 

Qualitative  Comments  on  the  Results 

The  numerical  method  presented  here  has  been  shown  to  be  a  very  efficient  and  com¬ 
petitive  technique  for  the  solution  of  the  3-D  Navier-Stokes  equations  over  hypersonic 
bodies.  The  iteration  count  for  the  solution  of  steady-state  problems  is  on  the  order  of 
one  hundred  for  the  test  cases  studied  here.  Further  work  needs  to  be  done  to  vectorize 
the  code  and  thus  reduce  the  CPU  time  per  iteration  per  grid  point  from  the  present 
9.8  x  10~4sec  on  a  Cray  X-MP  48. 

The  test  cases  have  demonstrated  several  important  features  of  the  implicit  Gauss- 
Seidel  line-relaxation  technique.  The  iteration  count  is  about  the  same  for  each  section 
of  the  X24C  test  case  as  for  the  biconic.  An  increase  in  the  number  of  points  in  the 
meridional  (f )  direction  does  not  increase  the  number  of  iterations  required  for  a  steady- 
state  solution.  The  mesh  spacing  near  the  wall  may  be  reduced  without  increasing  the 
number  of  iterations  or  the  CPU  time  needed  for  convergence. 

The  number  of  iterations  required  is  proportional  to  the  physical  length  of  the  sec¬ 
tion,  the  number  of  grid  points  in  the  i  (streamwise)  direction,  the  obliqueness  of  the 
shock  wave,  and  the  quality  of  the  guess  of  the  initial  bow  shock  wave  position.  The 
convergence  to  steady-state  is  somewhat  hindered  by  the  coupling  of  the  numerical 
method  to  the  present  shock-fitting  technique.  The  steady-state  solution  is  dependent 
on  having  the  correct  bow  shock  position  and  a  converged  flow  field,  thus  both  aspects 
of  the  solution  must  progress  in  time  at  approximately  the  same  rate  for  an  optimal 
convergence  rate  to  be  achieved.  It  appears  from  the  nature  of  the  solution  develop¬ 
ment  that  the  shock-fitting  routine  lags  the  solution  of  the  flow  field,  slows  down  the 
total  convergence  rate,  and  increases  the  number  of  iterations  required  for  a  steady- 
state  solution.  In  particular,  the  number  of  iterations  to  reach  a  steady-state  for  the 
biconic  at  a  =  0.0°  was  reduced  by  a  factor  of  two  by  specifying  the  converged  grid 
location  as  the  initial  guess  of  the  grid  shape.  This  suggests  the  need  for  an  improved 
shock-fitting  technique  or  a  different  approach  to  capturing  the  bow  shock  wave  to 


achieve  the  best  convergence  rate  for  this  numerical  method.  However,  in  spite  of  this 


inefficiency,  the  two  procedures  work  well  together  and  demonstrate  a  small  iteration 
count  for  a  steady-state  solution. 

Conclusions 

In  conclusion,  a  numerical  method  to  solve  the  three-dimensional  Navier-Stokes 
Equations  for  hypersonic  bodies  has  been  developed  and  tested.  The  test  cases  demon¬ 
strate  that  the  method  is  robust  and  gives  good  results  for  a  difficult  problem.  The 
technique  has  been  shown  to  demonstrate  a  reduction  in  the  number  of  time  steps 
required  to  achieve  a  steady-state  solution  of  approximately  two  orders  of  magnitude 
over  an  explicit  method.  The  primary  advantage  of  this  method  is  that  there  is  no 
penalty  for  making  the  mesh  finer  near  the  wall  to  resolve  viscous  and  thermal  effects. 
This  algorithm  represents  a  step  toward  being  able  to  compute  three-dimensional  flow 
fields  in  a  small  number  of  iterations  and  with  much  less  computer  time. 
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Figure  3.  Comparison  of  Surface  Pressure  Distribution  on  Biconic.  Experimental 
results  are  from  Coutler  and  Carver*,  a  =  0.0°,  M <*,  =  7.97,  Re  =  1.2  x  108/m. 
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Figure  5.  Temperature  Profiles  on  Biconic  Afterbody, 
a  =  0.0°,  =  7.97,  Re  =  1.2x  10d/m. 


Figure  7.  Surface  Pressure  Distribution  on  Biconic  Afterbody. 

a  -  7.0°,  \tao  -  7.07,  He  -  1.2*  108/rn. 


Figure  8.  X24C  Body  Surface 
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Figure  9.  Surface  Pressure  Distribution  on  Plane  of  Symmetry  for  X24C.  Experi¬ 
mental  results  from  Wannernwetsch9.  q  =  6.0°,  A/oo  =  5.95,  Re  =  16.4xl06/m. 


Surface  Pressure  Distribution  on  X24C  at  x/rn  =  108.  Experimental  results 
from  Wannernwetsch9.  a  =  6.0°,  M =  5.95,  Re  =  16.4  x  106/m. 


Transfer  Distribution  on  X24C  at  z/r„  =  108.  Experimental 
Vannernwetseh0  ami  Neumann1 1 .  n  -  6.0°,  A/oo  =  5.95.  Re  =  16. 4x 


Figure  1G.  Streamline  Traces  on  Leeward  Surface  of  X24C  at  x :r 
a  =  G  O  .  .\/r< j  -  5.95,  Re  =  IG.-4  x  10"  rn. 
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